Interplay of superconductivity and spin density wave order in doped graphene 
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We study the interplay between superconductivity and spin density wave order in graphene doped 
to 3/8 or 5/8 filling (a Van Hove doping). At this doping level, the system is known to exhibit weak 
coupling instabilities to both chiral d + id superconductivity and to a uniaxial spin density wave. 
Right at van Hove doping, the superconducting instability is strongest and emerges at the highest 
T c , but slightly away from van-Hove doping a spin-density-wave likely emerges first. We investigate 
whether at some lower temperature superconductivity and spin-density-waves co-exist. We derive 
the Landau-Ginzburg functional describing interplay of the two order parameters. Our calculations 
show that superconductivity and spin density wave order do not co-exist and are separated by 
first-order transitions, either as a function of doping or as a function of T. 



Two dimensional electron systems provide an ideal 
environment for exploration of many body physics. 
Graphene, as a new two dimensional electron system, 
may allow us to access new many body phases that have 
not been hitherto observed. Unfortunately, undoped sin- 
gle layer graphene seems to be well described by a non- 
interacting model pQ , with the vanishing density of states 
suppressing interaction effects. In order to access many 
body physics in graphene, one must sidestep the vanish- 
ing density of states. One way to do this is by doping. 
When graphene is doped to the M point of the Bril- 
louin zone - a doping level that corresponds to 3/8 (or 
5/8) filling (undoped graphene corresponds to 1/2 fill- 
ing) - the Fermi surface undergoes a topological transi- 
tion from a two piece to a one piece Fermi surface [2]. 
Associated with this topological transition is a divergent 
density of states, which gives rise to weak coupling insta- 
bilities to unusual many body states. Doped graphene 
thus provides a promising playground for exploration of 
new quantum many body states. 

The recent success of experimental efforts to dope 
graphene to the M point [3] has inspired a flurry of 
theoretical works studying many body physics in doped 
graphene jlHE] ■ It has been established [2] that the prin- 
cipal weak coupling instabilities are to chiral d+id super- 
conductivity and to a uniaxial spin density wave (SDW), 
with the superconducting instability leading at Van Hove 
doping [4 , and the SDW leading somewhat away from 
Van Hove doping |llj . It is not known, however, whether 
these two orders are mutually exclusive, or whether they 
co-exist in some range of temperatures and dopings. 

In this letter we demonstrate that for graphene near 
the M point, superconductivity and SDW magnetism are 
mutually exclusive orders. We derive Landau-Ginzburg 
action for two order parameters and show that the inter- 
play between quartic terms is such that the minimum of 
the action is when only one order parameter is nonzero. 
This result stands in stark contrast to pnictide materials, 
where Landau-Ginzburg analysis shows that supercon- 
ducting and SDW orders do co-exist [HI [16] . In doped 
graphene one expects to observe pure chiral superconduc- 



tivity at Van Hove doping, with a first-order transition to 
pure spin density wave order upon doping away from the 
Van Hove point. Our conclusions apply also to doped 
triangular lattice systems [T7], which have an identical 
low energy description near 3/4 filling. 

The model: Our point of departure is the tight binding 
model [2 , with the nearest- neighbor dispersion 



e k = ±tyi+4cos^cos 3 ^+4cos 2fc ^- M (1) 

where the overall sign is + or — depending on whether we 
are above or below half filling. For definiteness we take 
a plus sign. Van Hove doping then corresponds to fi = t, 
at which point the Fermi surface has the form shown in 
(Fig. [I]). The Fermi velocity vanishes near the hexagon 
corners Mi = (2tt/3,0), M 2 = (7r/3,7r/V3), M 3 = 
(— 7r/3, 7r/v3), which are saddle points of the dispersion: 

it 3t 

EfcRiMj = ^(3£ 2 -fc 2 ), £fcf«Af 2 = —2ky{ky - V3k X ), 

it 

£fcf«M3 = -j2ky(k y + V3k x ) (2) 

Each time k is a deviation from a saddle point. Saddle 
points give rise to a logarithmic singularity in the DOS 
and control physics at weak coupling. There are three in- 
equivalent nesting vectors Q a b connecting in-equivalent 
pairs of saddle points M a and M{, (see Figjl]): 

Qi = Q23 = ( 7r > 7r /v / 3); Q2 = Q31 = (tt, -tt/\/3) 
Q 3 = Q12 = (0,2tt/V3). (3) 

Each Qi is physically the same as — Qi because Qi is half 
of a reciprocal lattice vector. 

There are four different interactions between fermions 
near saddle points, i = 1 — 4, with momentum trans- 
fer near zero and near Qi. (Rcfs. HI EHH1 El- For our 
purposes, relevant interactions are density-density inter- 
action within one patch (54) and between patches (#2) 
and the interaction which describes hopping of a pair of 
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fermions from one patch to the other (33). The fourth in- 
teraction gi is the exchange interaction between patches. 
Interactions 32 and 33 renormalize particle-hole vertices 
and control the SDW instability, while interactions 33 
and 34 renormalize particle-particle vertices and control 
the superconducting instability (note that 33 contributes 



to both instabilities). 

The partition function Z of the model can be written 
as a functional integral over Grassman valued (fermionic) 
fields ip. We have Z = J D[ip,ip] exp ( — Sfy, ip)) , where 

S = J^ T £(k, t) (T is the temperature) and 



4>a,a(d T +E- k -H-g4lpa,ai'a,a)lpa, 



b=£a 




FIG. 1: The Fermi surface at Van Hove doping is a perfect 
hexagon inscribed within the hexagonal Brillouin zone. The 
hexagon has three inequivalent corners, labeled ¥1,2,3, which 
are saddle points of the dispersion and give rise to a divergent 
density of states. Each saddle point is perfectly nested with 
each other saddle point. The perfect nesting of the FS is 
broken only by third and higher neighbor hoppings, which 
are generally quite small. Meanwhile, the existence of saddle 
points is fully robust, being a consequence of a topological 
transition from a Fermi surface with two inequivalent pieces 
to a one piece Fermi surface 



Here a,b — 1, 2, 3 label which saddle point we are closest 
to, a and j3 are spin labels, and a is the opposite spin 
state to a. We have retained only those states that are 
close to the saddle points - this 'patch model' is exact in 
the limit of weak coupling [4] . 

This action displays instabilities towards d-wave 
superconductivity and SDW. We therefore decouple 
the interactions in the d-wave superconducting and 
SDW channels simultaneously, by means of two Hub- 
bard Stratanovich transformations. We introduce the 
Hubbard-Stratanovich superconducting fields A a = (33 — 
9i){'4'a,-\4 l a,i) ■ Since the superconductivity is known to 
be d + id [J, we set (Ai, A 2 , A 3 ) = A(l, e 2 "/ 3 , e^ 2 "/ 3 ) 
and describe superconducting fields by a single complex 
order parameter A. We also introduce the three SDW 
order parameters M a b = (32 + gaXVvtV'fc^)- Since the 
SDW order is known to be uniaxial we can re- 



place the three vector order parameters M12, M23, M31 
by a single scalar SDW order parameter M, which rep- 
resents the magnetic order along the SDW axis. Since 
the system has 0(3) spin rotation symmetry, the SDW 
axis can be chosen to coincide with the z axis with- 
out loss of generality. Finally, we introduce the Nambu 
spinor \a, a four component spinor defined according 
to Xa = (VVh^a4'^a,4.' _ VV|-)- The action after Hub- 
bard Stratanovich transformation can be written in the 
Nambu spinor basis as 



C = 



XaQ a b Xb 



Qaa 



'a^b 



M |A| 2 

.92+33 33-34 

(d T \ 2 + f>k - ju)er 3 + Acr + + A* a. 



(5) 



(6) 



Here the ai are Pauli matrices acting in the particle-hole 
space, the rji are Pauli matrices acting in the spin space, 
I2 is a two dimensional identity matrix, and a± = a\±a2- 
The notation we have used is borrowed from 18J. 

We can now integrate out the fermions exactly to ob- 
tain an action purely in terms of the superconducting 
and SDW order parameter fields, 



C = 



M 2 



92 +33 33 - 34 



-Trln0- 1 (A,M) (7) 



where the trace goes over Nambu spinor indices, and also 
over imaginary time and over momentum. We now define 
G to be the "bare" (matrix) Green function evaluated at 
A = 0, M = 0, and define matrix order parameters A 
and M, such that Q^ 1 = G 1 + A + M. We can then 
write 



Trln G- 1 (1 + G(A + M N 



constant + Tr In 1 + G(A + M) . (8) 



Tr ln£" 



It is convenient to explicitly write out the expressions 
for G, A and M. We adopt the shorthand ^(u^k) = 
1 / (iui n ± (ek — fi)j, where u) n = (2n + 1)ttT are fermionic 



3 



Matsubara frequencies. Using the shorthand, we can de- fine the various matrices as 
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(9) 



Thus far, everything we have done has been exact. We now work close to T c and perform a double expansion of ([8]) 
in small |A| and small M. We terminate the expansion at quartic order in both fields and drop all terms that are 
odd in powers of M or A as they vanish upon taking the trace. We then obtain for the Tr In term the expression 



Tr 



- (GAGA + GMGM) - (GAG AGAGA + GMGMGMGM + 4GMGMAG A + 2GMGAGMGA) 



We have made use of the fact that the trace of a product of matrices is invariant under a cyclic permutation of the 
matrices. Evaluating the traces and substituting back into ([7| leads to the expression 



£ = ai(T - T C )|A| 2 + a 2 (T - T N )M 2 
where we have defined the expansion coefficients 

d 2 k 



KAAr + K 2 M + 2K 3 \A\ 2 M 2 



K 2 =3T^Jj§^ 



(F+K,k)F-K,k)f 
F+K,k) 2 F+K,k + Q 1 ) 2 



K 3 =6T£ Wn f I g j ^ F+K,k)F-K,k) J F+K,k + Q 1 ) J F-K,k + Q 1 )cos 

k) 2 F-(w„, k)F+K, k + QO + (F+ o F~) 



2F+K,k) 2 F+K,k + Q 1 )F+K,k- 

Y 



Q 2 ) + (^ 



(10) 

(ii) 



r 



Terminating the expansion at quartic order in both or- 
der parameters is justified if the quadratic terms for su- 
perconducting and SDW order change sign at about the 
same critical temperature, T c T/y. Renormalization 
group analysis shows that the couplings which determine 
both superconducting and SDW instabilities diverge at 
the onset of the first instability upon lowering T [?J [71 [3] • 
The critical temperatures, however, are determined by 
superconducting and SDW susceptibilities, which gener- 
ally have different exponents (different anomalous dimen- 
sions). It has been demonstrated on general grounds [19] 
that different orders emerge simultaneously when their 
anomalous dimensions r\i > 1. In our case, the anomalous 
exponents for superconducting and SDW susceptibilities 
have been calculated in [2] . Using results from that work, 
we find that r]sc — 1-48 > 1, while tjsdw = 0.97 (the 



results are for perfect nesting). Because tjsdw < 1> 
T N < T c . However, since ijsdw is ver Y close to one, 
we expect that Tn in (10) is only slightly lower than T c , 



in which case the expansion up to quartic order in A 
and M is justified. For dopings slightly away from the 
van-Hove one, we expect SDW to be the leading instabil- 
ity [7J H3 HU [T2] . In this case, T N > T c , again we assume 
that the difference between the two critical temperatures 
is small. 

It was shown in the context of pnictides [13 [H] that a 
free energy of the form ( 10 1 leads to co-existence of the 
two order parameters if K\K 2 — K 2 > 0. We computed 
the coefficients K\, K 2 and K 3 in our case by explicitly 
integrating over fermionic momenta and summing over 



fermionic frequencies in ( 12 ) (for details see Supplemen 



4 




FIG. 2: Diagrammatic representation of the two processes 
that couple superconductivity and magnetism at quartic or- 
der in the Landau- Ginsburg expansion. Lines with arrows 
represent fermion propagators, dotted lines represent SDW 
order parameter M, dashed lines represent superconducting 
order parameter A. 



tary material). We obtain, with logarithmic accuracy 



L05 , * ,. 
— = In — + subheading 



!- 05 , t ,, 

K 2 = —7= In — + subleadmg = K\ 

1.05 ( cos(^) + 2 
Ka = — In — + subleading 

= i^cos(^) + 2). (12) 

Note that there are two processes which contribute to the 
coefficient K 3 . The processes are represented diagram- 
matically in Fig.Q. The process shown in Fig.|2]a) is 
sensitive to the chirality of the superconducting order 
parameter, because of the dependence on the phase dif- 
ference between different saddle points, and gives rise 
to the cos(^) term. Because cos(27r/3) = —1/2, this 
process gives rise to an effective attraction between su- 
perconductivity and spin density waves. This effective 
attraction is, however, outweighed by a larger (chirality 
independent) repulsion between the two order parame- 
ters, coming from the processes shown in Fig.([2jb). The 
prefactors in our case are such that K 3 = \K X > 0. Com- 
paring K\Ki and iff we see that in the case of doped 
graphene K\K 2 — K\ < 0, so that co-existence is dis- 
favored. The system only allows one order parameter 
to exist, even when T c = Tjy. A direct second order 
transition between superconducting and SDW orders is 
Landau forbidden, since the symmetry group of one or- 
dered phase is not a subgroup of the symmetry group of 
the other ordered phase, and as a result the transition 
separating the region when A / 0,M = from the re- 
gion where A = Q, M ^ is expected to be first order. 
( Although we cannot exclude a non-Landau continuous 
transition between the two ordered states.) 

The fact that T c ^ Tn makes co-existence even less 
likely. We therefore conclude that there is no co-existence 
of superconducting and SDW order in doped graphene. 

In pnictides, the structure of Ki,K 2 , K 3 is quite simi- 
lar [TB] (modulo that there is no log t/T c term), but the 
argument of cos in K 3 is the phase difference between the 



gaps on hole and electron FSs. For superconductiv- 
ity, the argument is tt, in which case K 3 — K\ = K 2 . 
Then, K\K% = K 3 , and one has to include sublead- 
ing terms to verify whether the two orders can co-exist. 
The subleading terms are the ones which break the nest- 
ing between hole and electron pockets, and the analysis 
shows [13 US] that s H superconducting and SDW or- 
ders co-exist in some range of parameters. In graphene, 
the argument of cos is 27r/3, and such co-existence does 
not occur. 

To conclude, we have demonstrated that superconduc- 
tivity and spin density wave order are mutually exclusive 
in graphene doped near the M point of the Brillouin zone 
(a Van Hove doping) . Sufficiently close to the Van Hove 
point, we expect to see pure chiral superconductivity, and 
somewhat away from the Van Hove point we expect to see 
pure spin density wave order. The results stand in stark 
contrast to pnictides, where there can be co-existence be- 
tween spin density waves and superconductivity. 

We acknowledge useful conversations with C. Batista, 
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J. Schmalian, and R. Thomale. A.V.C. is supported by 
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APPENDIX 



Evaluating K\ 

We start with the expression 



Where the sum goes over all Matsubara frequencies, and the integral goes over all wavevectors k up to a UV cutoff of 
|k| ~ 0(1), at which point the dispersion relation changes. The UV cutoff must be retained in the integrals, because 
the integrals are log divergent in the UV if we ignore the cutoff (the integrals are convergent in the infrared at nonzero 
temperature). We now scale out uj n , and define the rescaled co-ordinates x = y3i/ (4bj n )k x , V — y/3t/(4u} n )k y . The 
expression for K\ can then be recast as 



/t/T a 

z - 3irH\u„\ 3 L^JTJt c _ (l + (3x 2 - y 2 ) 2 )' 



We now integrate first over — oo < y < oo, and then over —^JtjT c < x < y/t/T c , and expand the resulting expression 
to leading order in large t/T c (the manipulations are all done on Mathematica) . We obtain the result 

Kl = Tc ^ ^RF ( 27! ln k + subleadin s) ■ ( 15 ) 

The same result is obtained, with logarithmic accuracy, if we first integrate over — oo < x < oo, and then over 
-\/t/T c < y < We now recall that ui n = {In + l)vrT c , and that |n+ ^ /2|3 = 14^(3) w 16.8 (the sum may 

again be evaluated on Mathematica). Thus we obtain, with logarithmic accuracy, the result quoted in the main text, 
i.e. 

14C(3) , t 1.05 t 

K 1 = — ln — w -= ln — (16) 

16V3vr 4 tT 2 T c y/SwHT^ T c v ; 



Evaluating K2 

The co-efficient K 2 was evaluated already in [TT]. We can write K 2 = 6Z1 + 12Z 2 , where Zi and Z 2 are coefficients 
that were defined in [TT] and calculated in the supplement to TJJ. (Note that there is an overall factor of 2 relative 
to |11| , which comes about because we have doubled the number of degrees of freedom in going to the Nambu spinor 
representation. However, this overall factor of 2 multiplies all terms in our free energy, and thus has no physical 
significance). In [TT], it was shown that the term Z\ was larger than Z 2 by a factor of ln t/T c , which is a large 
number at weak coupling. Thus, we can neglect Z 2 with logarithmic accuracy, and say K 2 = 6Z1. The co-efficient Z\ 
was calculated in |11) . however, the calculation there had a factor of 2 error, which was unimportant for the physics 
considered in [TT] but is important here. Therefore, we redo the calculation of Z\. 

We wish to evaluate 

Z i= T z2j ^^ + (k, W „) 2 F+(k + Q, Wn ) 2 (17) 

The integral over the Brillouin zone is dominated by those values of k where both Green functions correspond to 
states near a saddle point. Expanding the energy about the saddle points, we rewrite the integral as 



Zi^rr/^S 5-^ 9 (18) 

^ ' (2^) 2 _ 3|a ( 3 fe2 _ fc 2 )) 2 ^„ - ^2k y (k y - V3k x )) 2 



d 2 k 1 



() 



Where the integral is understood to have a UV cutoff for k of order 1. We now define a = \f3t\/A(k y — \/3k x ) and 
b = y/3ti/4(k y + V3k x ), and rewrite the above integral as 

t"' 3V3!i J-vt; W 



r ) 2 (iw„ + abj (iu> n — a(a + 6))" 



We now define x = ab and rewrite the integral as 



z =T Y^- f -- I - 1 (20) 

We now assume Tjy <C t\ (which should certainly be the case for weak/moderate coupling). In this limit, we can 
perform the integral over x approximately, using the Cauchy integral formula, to get 

= T y 2 da 1 2»signa;» =T y" 4 Z' 7 * 1 ^ 1 »signcj n (q 2 + 2icj») 3 

1 tr3V3t 1 1^2^|a|(a2- 2lW „) 3 tTS^J-Vtr^W (« 4 + 1 j 

The imaginary part of the above integral is odd in uj and hence vanishes upon performing the Matsubara sum to leave 
an integral that is purely real 



z =T sp 8 KI da 1 4^ -3a 4 ^ 8\uj n \ daj__ 

f - 3^3*1 J-VE 2tt \a\ (a 4 + 4^)3 ~ 3^ ^ 27r | a | (a 



(22) 



' ' x 4 W 2)3 

with logarithmic accuracy. Performing the integral over a (again with logarithmic accuracy) gives 

El 1 U 1 , U , 16.8 In 4± 
7^1 ^ln— = (16.8 In— + 10.5) « (23) 
__. n 127rV3ti Kl 3 w„ Q&n*y/VI%tx V 2ttT ; 96^^^^ V ; 

Where we take ui n = 2n(n + l/2)T/v, T — T/v and perform the discrete sum on mathematica. The error in the 
supplement to lT| was in the last line of the calculation. 



Evaluating K3 

There are two distinct contributions to K%, and we evaluate both in turn. We can write K 3 = K% + K^, where 

Kl =6T-£ Un J (0 F+K,k)F-K,k)F+K,k + Q 1 )F-K,k + Q 1 )cos(0 k -^ k+Q ) 

K b 3 =6TE w „/(0 ^ + K,k) 2 F'K,k)^ + ( W „,k + Q 1 ) + (F+oF-) (24) 

The first contribution K 3 comes from processes of the form shown in Fig.([2ja), and is sensitive to the chirality of the 
superconducting order parameter (it depends on the difference in the phase of the superconducting order parameter at 
different points on the Fermi surface). This process leads to an attraction between chiral superconductivity and spin 
density waves The second contribution K\ comes from processes of the form shown in Fig.j2]b), and is insensitive 
to the chirality of the superconducting order parameter. This process leads to a repulsion between any kind of 
superconductivity and spin density waves. The second process dominates (because of purely numerical prefactors), 
so superconductivity and spin density waves do repel - but the repulsion is too weak to prevent co-existence. 
Let us first calculate iff. For d + id pairing, we have 6>k — $k+Q = 47r/3. Thus, we have 



KS =6T c £_J" 7 f^ F+(u n ,k)F-(u n ,k)F + (u n ,k + Qi)F-{u n ,li + Qi)coa 



4?r 

(2ir) 2 x ' y^UT^-J-L \w n , MX.JJ.- \<M n , xv ~r v^i \u/ n) r>. -p v^jj wo — 

= 6T c cosf X u Jt$& ^ ~ ^ (3fc| _ k2y)2) ^ " ^ Ak2y{ky _ ^ kx)2) ■ (25) 



We scale out us n and define rescaled variables x — y // 3t/(4 : uj n )k x , y = ^3t/ {Auj n )ky. The expression for K% can then 
be recast as 

4tt\^ 2T c fV*/T° dxdy 



A;; 3 J f-7r2t| w „| 3 7.^(1 + (3x2- y2)2) (l + V^-^) 2 )' ( - (,) 
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We define the new coordinates a = y — \/3x, b = y + y/3x, and hence re-express the above integral (with logarithmic 
accuracy) as 

K a _ ( 4n \sr Tc dadb 

~ V V tT V3n*t\u n \* J_ VWc (l + a 2 & 2 )(l + a 2 (a + fo) 2 )' (2?) 

We integrate over —y/t/T c < b < \Jt/T c (on Mathematical and expand the resulting expression to leading order in 
large t/T c . This leads to the expression 



K3 = cos ( ^ ) y — ^ ( [ da—^—j- ) (28) 



It should be remembered that the expansion in large t/T c is valid only for a 2 t/T c 1, thus the above integral 
implicitly carries an infrared cutoff on the scale a w \JtjT c . Performing the integral with this infrared cutoff, we 
obtain the expression 

Performing the summation over n on Mathematica, as before, we obtain 

14C3 . t -1 16.8 t -1.05 t 1 

AT" = cos — ^ In — rts = In — = — = In — = --K\ 

3 V 3 / mVSnHT? T c 2 16\/37r 4 iT 2 T c <2^HT? T c 2 



(30) 



Note the crucial minus sign that comes from the chirality sensitive cos factor - this particular term represents an 
attraction between magnetism and chiral superconductivity. 
We now turn our attention to the second term, K%. We have 

K\ =6T c E w „/(0 F+K,k) 2 F"K,k)F+(a;„,k + Q 1 ) + (F+ F~) 

(iuj n + f (3fc 2 - fc 2 )) (iw„ + ^2k y (ky - V3k x )) + (w„ -> -w„) 



d"fc 

( 27r ) 2 A .2 , 9* 2 /o/,2 7„2\2\ 2 / 



(a; 2 + f (3fc 2 - fc 2 ) 2 )> 2 + % 4k*{ky V3k x y) 

= 12T V f ^ fc ^ ~ ^ 3fe ^ ~ k y^ 2k y( k v ~ / 31 x 

( " )2 K + 5f(3fc 2 -fc 2 ) 2 ) 2 ( W 2 + 5f4fc 2 (fc y -^fc x ) 2 ) 



Again, we scale out w„ and define the rescaled variables x = \/3t/ (4u n )k x , y — ^/3t/(4w n )fc 3/ . The expression for K\ 
can then be recast as 

^ = E -T^ll ^ l-(^- y p2y(y-V3x) 

^ 7r 2 i|o; Tl | 3 J^^/TJ¥ C (l + (3x 2 -y 2 ) 2 ) (l + 4y 2 (y-\/3x) 2 ) 

We define the new coordinates a = y — V~3x, b = y + \/3x, and hence re-express the above integral (with logarithmic 
accuracy) as 

Kt = V r 2T < dadb ^±^±3 (33 ) 

^ y/torH\w n \* J-^JT C (l + a 2 6 2 ) 2 (l + a 2 ( a + 6) 2 ) 

We integrate over —y/t/T c < b < \ft/T c (on Mathematica), and expand the resulting expression to leading order in 
large t/T c . This leads to the expression 

^ = ?^bF(/ da 2|a|(l + aV4)0 (34) 
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It should be remembered that the expansion in large t/T c is valid only for a 2 t/T c 3> 1, thus the above integral 
implicitly carries an infrared cutoff on the scale a w \JtjT c . Performing the integral with this infrared cutoff, we 
obtain the expression 

Kl = y_^L ln-= ? m In - = 2^ (35) 

Putting things together, we have 

K 3 = K$ + Kl = ^(cos y +2) = + 2) = ^! (36) 

quoted in the main text. 



